#!/usr/bin/python
"""
This program pick primers for a list of given regions.

Usage: dissertation_PrimerPicker.py genome.fasta regions.list
Note: Need to run in primer3 program folder:
/home/js/tools/primer3-2.3.2

Regions need to be in tab-delimited format like this:

region1	1000	2000
region2	6000	10000
region3	12000	14000

Currently, the values are hard-coded for long PCR primers to produce
product size of about 15kb.
Should check the results manually to pick best primers but eventually this script
should be able to pick the best one for given conditions (TODO).

NOT DONE YET. NEED TO CONTINUE.

"""


import os
import sys
import re
from Bio import SeqIO
from Bio import Motif
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq

primers = re.compile('   \d+ (\w+)\s+\d+ \d+\s+ 0 \d+.*')

#from subprocess import call

#Def to write primer3 config files
def MakePrimerConfig(seqid, seq, left, prodsize):
    parameters = []
    parameters.append("SEQUENCE_ID=" + seqid)
    parameters.append("SEQUENCE_TEMPLATE=" + seq)
    parameters.append("SEQUENCE_TARGET=" + str(left) + "," + str(prodsize))
    parameters.append("PRIMER_TASK=pick_detection_primers")
    parameters.append("PRIMER_PICK_LEFT_PRIMER=1")
    parameters.append("PRIMER_PICK_INTERNAL_OLIGO=0")
    parameters.append("PRIMER_PICK_RIGHT_PRIMER=1")
    parameters.append("PRIMER_OPT_SIZE=25")
    parameters.append("PRIMER_MIN_SIZE=24")
    parameters.append("PRIMER_MAX_SIZE=28")
    parameters.append("PRIMER_MIN_GC=50")
    parameters.append("PRIMER_MAX_GC=65")
    parameters.append("PRIMER_MAX_END_GC=2")
    parameters.append("PRIMER_MAX_NS_ACCEPTED=0")
    parameters.append("PRIMER_PRODUCT_SIZE_RANGE=14500-15500")
    parameters.append("PRIMER_THERMODYNAMIC_PARAMETERS_PATH=./src/primer3_config/")
    parameters.append("P3_FILE_FLAG=1")
    parameters.append("PRIMER_EXPLAIN_FLAG=0")
    parameters.append("=")
    outhandle = seqid + ".txt"
    outfile = open(outhandle, "w")
    outfile.write("\n".join(parameters))

genome = SeqIO.read(sys.argv[1], "fasta")

regionsfile = open(sys.argv[2], "rU")
regions = regionsfile.readlines()

#parsing the regions to pick primers from
for line in regions:
    l = line.split('\t')
    name = l[0]
    start = int(l[1])
    stop = int(l[2])
    leftflank = 400
    productsize = 14500
    #productsize = (stop - start + 1)
    sequence = str(genome.seq[start-100:stop+100])
    MakePrimerConfig(name, sequence, leftflank, productsize)
    primerout = name + ".txt"
#    call(["primer3_core", "<", primerout])
    os.system("primer3_core < " + primerout)
    forprimer = open(name + ".for", "r")
    forward = forprimer.readlines()
    revprimer = open(name + ".rev", "r")
    reverse = revprimer.readlines()
    if len(forward) > 3:
        for i in forward[4:]:
            mtf = Motif.Motif(alphabet=IUPAC.unambiguous_dna)
            tmp = i.split(' ')
            pseq = tmp[4]
            #NOT DONE YET. NEED TO CONTINUE
    else:
        print "Forward primer not found. Edit the config file to manually design."

    if len(reverse) > 3:
        for j in reverse[4:]:
            mtf = Motif.Motif(alphabet=IUPAC.unambiguous_dna)
            pattern = 


#Sifting through the results to pick the best pairs
for line in regions:


regionsfile.close()
